// pt16.C -- post ta pid going on -- z = 16

// p[0]: aoz; x[0]: beta
// return p/Z = poz(beta, aoz);
double pv(double *x, double *p);
double pbg(double *x, double *p);

void pt16(){
  TFile *f = new TFile(getenv("TMP1"), "UPDATE");
  TTree *etf = (TTree *)f->Get("etf");
  // auto f = gDirectory;
  // TChain *etf = new TChain("etf");
  // etf->Add("kr0.root");
  // etf->Add("ch.root");
  // etf->Add("c.root");
  cout << "Totally " << etf->GetEntriesFast() << " events" << endl;
  TCutG *az_[3] = {(TCutG *)f->Get("s33"), (TCutG *)f->Get("s34"), (TCutG *)f->Get("s35")};
  TCut cut = "1 "
    // "&& x_f1>=-17.&&x_f1<=-11. "
    // "&& pdcTaUWire[1][0]>=-15.&&pdcTaUWire[1][0]<=-11. "
    // "&& pid.betaz[0]()>0.6&&pid.betaz[0]()<0.7 "
    // "&& fabs(globalTracker.dx2())<4. " // good pid
    "&& pid.poz()>1300. "
    "&& sqrt((pTAU[0]+7.)*(pTAU[0]+7.)+(pTAU[1]+1.5)*(pTAU[1]+1.5))<25. " // within TA
    "&& fabs(pT1[0]+6.2)<24. && fabs(pT1[1]+3.2)<24. " // within tstop
    "&& z0c>=15.5&&z0c<=16.5 "
    "&& (s33 || s34 || s35) "
    // "&& tof1>=135.&&tof1<=155. "
  ;
  static const double zr21[6] = {1.9, 1.97, 2.018, 2.069, 2.1186, 2.159}; // z=21
  static const double zr20[6] = {1.9, 1.97, 2.02, 2.075, 2.123, 2.2}; // z=20
  static const double zr19[6] = {1.9, 1.974, 2.023, 2.078, 2.133, 2.2}; // z=19
  static const double zr18[4] = {2.024, 2.077, 2.134, 2.193}; // z=18
  static const double zr17[6] = {1.9, 1.97, 2.027, 2.0845, 2.15, 2.216}; // z=17
  static const double zr16[6] = {1.9, 1.97, 2.027, 2.091, 2.15, 2.214}; // z=16
  char cutaz[128]{};
  const double zz = 18.;
  // auto zr = zr18;


  const char *az[3] = {"s33", "s34", "s35"};
  // TCut az[3] = {"s33", "s34", "s35"};
  // for(int i = 0; i < 3; i++){
  //   sprintf(cutaz, "aoz0c>%1.3f && aoz0c<%1.3f", zr[i], zr[i+1]);
  //   az[i] = cutaz;
  // } // end for over az cuts
  // sprintf(cutaz, "aoz0c>%f", zr[5]); az[5] = cutaz;
  etf->Draw(">>el", cut, "entrylist");
  TEntryList *el = (TEntryList *)f->Get("el");

  etf->Draw("z0c:aoz0c>>h1(1000, 1.9, 2.4, 1000, 14., 39.)");
  etf->SetEntryList(el);
  etf->Draw("et.tof2:aoz0c>>h2(100, 1.9, 2.4, 150, 60., 80.)");
  etf->Draw("pid.aoz()>>h3(103, 1.5, 2.6)");
  // etf->Draw("et.tof2:et.tof1>>h2(200, 138., 148., 200, 61., 74.)");
  // etf->Draw("pid.beta():pid.betaz[0]()>>h3(200, 0.59, 0.62, 200, 0.51, 0.64)");
  // etf->Draw("postMagTrackerX.GetHitPos(\"tofwB\"):aoz0c>>h4(200, 1.9, 2.4, 200, -700., 700.)");
  // etf->Draw("pid.aoz():aoz0c>>h4(100, 1.9, 2.4, 120, 1.9, 2.6)");
  etf->Draw("pid.bg0():pid.bg()>>h4(500, 0.5, 0.9, 500, 0.5, 0.9)", "aoz0c>0.", "col");
  // etf->Draw("postMagTrackerX.ang()-globalTracker.GetPostTaTrack().ang():aoz0c>>h4(200, 1.9, 2.4, 200, 10., 20.)");
  // etf->Draw("postMagTrackerX.GetHitPos(\"tofwB\")>>h6(300, -700., 700.)");
  etf->Draw("pid.poz():pid.beta()>>h5(100, 0.5, 0.9, 100, 1100., 1600.)", az[0]);
  etf->Draw("pid.aoz()>>h6(103, 1.5, 2.6)");

  TCanvas *c = new TCanvas("c", "Tstop Statistics", 1680, 1050);
  // TCanvas *c = new TCanvas("c", "Tstop Statistics", 1680, 1050);
  c->Divide(3,2); char h[64]; TH1 *b;
  for(int i = 0; i < 6; i++){
    sprintf(h, "h%d", i+1); f->GetObject(h,b);
    c->cd(i+1)->SetGrid();
    if(2 != i) b->Draw("col");
    else b->Draw("e");
  } // end for

  c->cd(1);
  etf->Draw("z0c:aoz0c>>h1_0(1000, 1.9, 2.4, 1000, 14., 39.)", "", "same");

  c->cd(6);
  etf->Draw("pid.aoz()>>ha0(103, 1.5, 2.6)", az[0], "same");
  etf->Draw("pid.aoz()>>ha1(103, 1.5, 2.6)", az[1], "same");
  etf->Draw("pid.aoz()>>ha2(103, 1.5, 2.6)", az[2], "same");
  // etf->Draw("pid.aoz1()>>ha3(103, 1.5, 2.6)", az[3], "same");
  // etf->Draw("pid.aoz1()>>ha4(103, 1.5, 2.6)", az[4], "same");
  double aoz[5]{};
  double par[9]{}; // gaussian for the partial aoz-s
  for(int i = 0; i < 3; i++){
    sprintf(h, "ha%d", i); f->GetObject(h, b);
    b->SetFillColor(i+1);
    b->Fit("gaus", "Q");
    TF1 *gaus = b->GetFunction("gaus");
    par[3*i] = gaus->GetParameter(0); // constant
    aoz[i] = par[3*i+1] = gaus->GetParameter(1); // mean
    par[3*i + 2] = gaus->GetParameter(2); // sigma
  } // end for over i
  // etf->Draw("pid.aoz()>>ha5(103, 1.9, 2.6)", az[5], "same");
  // ((TH1F *)f->Get("ha5"))->SetFillColor(5);
  c->cd(3);
  TF1 *aoz_f1[3];
  for(int i = 0; i < 3; i++){
    sprintf(h, "aoz_f1%d", i);
    aoz_f1[i] = new TF1(h, "gaus", 1.5, 2.6);
    // aoz_f1[i]->SetLineStyle(2);
    aoz_f1[i]->SetLineColor(i+1);
  } // end for over i
  TF1 *aoz_f1_merge = new TF1("aoz_f1_merge", "gaus(0)+gaus(3)+gaus(6)", 1.5, 2.6);
  aoz_f1_merge->SetLineColor(4); aoz_f1_merge->SetLineWidth(4);
  aoz_f1_merge->SetParameters(par);
  // for(int i = 3; i--;) aoz_f1_merge->FixParameter(3*i+2, par[3*i+2]);
  for(int i = 3; i--;) aoz_f1_merge->SetParLimits(3*i+2, 0.01, 0.03);
  ((TH1F *)f->Get("h3"))->Fit(aoz_f1_merge);
  aoz_f1_merge->GetParameters(par);
  for(int i = 0; i < 3; i++){
    aoz_f1[i]->SetParameters(&par[3*i]);
    aoz_f1[i]->Draw("same");
  } // end for over i


  c->cd(5);
  etf->Draw("pid.poz():pid.bg()>>hp0(100, 0.5, 0.9, 100, 1100., 1600.)", az[0]);
  etf->Draw("pid.poz():pid.bg()>>hp1(100, 0.5, 0.9, 100, 1100., 1600.)", az[1]);
  etf->Draw("pid.poz():pid.bg()>>hp2(100, 0.5, 0.9, 100, 1100., 1600.)", az[2]);
  // etf->Draw("pid.poz1():pid.beta1()>>hp3(100, 0.5, 0.9, 100, 1100., 1600.)", az[3]);
  // etf->Draw("pid.poz1():pid.beta1()>>hp4(100, 0.5, 0.9, 100, 1100., 1600.)", az[4]);
  // etf->Draw("pid.poz():pid.beta()>>hp5(100, 0.5, 0.9, 100, 1100., 1600.)", az[5]);
  for(int i = 0; i < 3; i++){
    sprintf(h, "hp%d", i); f->GetObject(h, b);
    b->SetMarkerStyle(3); b->SetMarkerColor(i+1);
    b->Draw("same");
  } // end for
  TF1 *faz[5];
  for(int i = 0; i < 3; i++){
    sprintf(h, "pv%d", i);
    // aoz[i] = (zz*2+1+i) / zz;
    faz[i] = new TF1(h, pbg, 0.4, 0.8, 1, 1);
    // faz[i] = new TF1(h, pv, 0.4, 0.8, 1, 1);
    faz[i]->SetParameter(0, aoz[i]);
    faz[i]->SetLineColor(i+1);
    faz[i]->Draw("same");
  } // end for

  c->Print("pt16.pdf");
  cout << "\033[36;1m" << cut << "\033[0m" << endl;
} // end pt

// p[0]: aoz; x[0]: beta
// return p/Z = poz(beta, aoz);
double pv(double *x, double *p){
  static const double cc = ETFMath::CC(), c0 = ETFMath::c0();
  const double aoz = p[0], beta = x[0];
  const double brho = aoz/cc*ETFMath::BetaGamma(beta);
  const double poz = brho * c0;
  return poz;
} // end function pv
double pbg(double *x, double *p){
  static const double cc = ETFMath::CC(), c0 = ETFMath::c0();
  const double aoz = p[0], bg = x[0];
  const double brho = aoz/cc*bg;
  const double poz = brho * c0;
  return poz;
} // end function pv
